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We have considered two classical lattice-gas models, consisting of particles that carry multicomponent mag- 
netic momenta, and associated with a two-dimensional square lattices; each site can host one particle at most, 
thus implicitly allowing for hard-core repulsion; the pair interaction, restricted to nearest neighbors, is ferro- 
magnetic and involves only two components. The case of zero chemical potential has been investigated by 
Grand-Canonical Monte Carlo simulations; the fluctuating occupation numbers now give rise to additional 
fluid-like observables in comparison with the usual saturated-lattice situation; these were investigated and their 
possible influence on the critical behaviour was discussed. Our results show that the present model supports a 
Berezinskil-Kosterlitz-Thouless phase transition with a transition temperature lower than that of the saturated 
lattice counterpart due to the presence of "vacancies"; comparisons were also made with similar models studied 
in the literature. 
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I. INTRODUCTION 

Planar rotators are models with potential interactions in- 
volving two-component spins; on the other hand, in the XY 
model, the spins have three components, only two of which 
are involved in the interaction, and the XY model can be re- 
garded as an extremely anisotropic (easy-plane) Heisenberg 
model. Actually, the terminological convention adopted here 
as well as by other Authors is not always followed, and the 
name "XY model" is sometimes used in the Literature to in- 
dicate planar rotators; on the other hand, both models are 
known to produce the same universality class. The two named 
models have been extensively studied, and possess a rich vari- 
ety of applications in Statistical as well as Condensed Matter 
Physicsi2i2. 

These models have been especially studied in their 
saturated-lattice (SL) version, where each site is occupied by 
a particle; as for symbols, classical SL spin models are first 
defined here: we consider a classical system, consisting of 
«-component unit vectors (usually n = 2,3), associated 
with a cZ-dimensional (bipartite) lattice Z d (mostly d = 2); let 
Xk denote dimensionless coordinates of the lattice sites, and 
let Uk, a denote Cartesian spin components with respect to an 
orthonormal basis e ff ; particle orientations are parameterized 
by usual polar angles {ifj} (« = 2) or spherical ones {(8j, <pj)} 
(» = 3); 

The interaction potentials considered here are assumed to 
be translationally invariant, restricted to nearest-neighbours, 
ferromagnetic (FM); they are, in general, anisotropic in spin 
space, and possess (9(2) symmetry at least: 

W jk = -e£l Jk , (1) 

where 

Q jk = au pjUk,,, + b^^u j_ a Uk y a', 

e > 0, < a < b, b > 0, max(a,Z?) = 1. (2) 



Here e is a positive quantity setting energy and temperature 
scales (i.e. T = kst/e, where t denotes the absolute tempera- 
ture); the corresponding (scaled) Hamiltonian reads 

A = - J] Cljk, (3) 

U<k) 

where }Z{j<k} i s restricted to nearest neighbours, with each dis- 
tinct pair being counted once. When d = 1,2, potential mod- 
els defined by Eq. (O produce orientational disorder at all 
finite temperatures (thus no orientational ordering transition 
at finite temperature)^; in the following, the discussion shall 
be essentially restricted to the two cases n — 2, a = b — 1 
(PR) and n = 3, a = 0, b = 1 (XY), respectively. 

It is by now well known that the PR model defined by 
d — 2 supports a transition to a low-temperature phase with 
slow (inverse-power-law) decay of the correlation function 
and infinite susceptibility; this is the very extensively studied 
Berezinskii-Kosterlitz-Thouless (BKT) transition (see, e.g. 
RefSi 1,2 ' 3 ' 6,7,8 and others quoted therein), whose existence was 
proven rigorously by Frohlich and Spencer—: the transition 
temperature was estimated to be ®pr = 0.88 + 0.01—, and a 
more recent and more refined result is ®pr - 0.8929^; the 
specific heat exhibits a maximum at a higher temperature^. 
Moreover, anisotropic models defined by d = 2, n — 3, < 
a < b have been studied as well and proven to support a BKT 
transition when the ratio a/b is sufficiently small 10 ; an esti- 
mate of the transition temperature for the XY case is given by 
®xy = 0.700 ± 0.005iii^. 

Lattice-gas (LG) extensions of the SL models can be de- 
fined as well, where each lattice site hosts one particle at most, 
and site occupation is also controlled by the chemical poten- 
tial //; such models have often been used in connection with 
alloys and absorption; this methodology somehow allows for 
pressure and density effects. Lattice-gas (LG) extensions of 
the continuous-spin potential model considered here are de- 
fined by Hamiltonians 

U<k) k 



2 



where Vk — 0, 1 denotes occupation numbers; notice that A < 
reinforces the orientation-dependent term, whereas A > op- 
poses it, and that a finite value of A only becomes immaterial 
in the SL limit ji — > +oo. 

For a square lattice, the SL-PR model produces a low- 
temperature BKT transition; the existence of such a transition 
for the LG counterpart has been proven rigorously as welU 4 .. 
In addition to the named rigorous results, still comparatively 
little is known about the above magnetic LG models with con- 
tinuous spins, in marked contrast to the vast amount of infor- 
mation available for their SL counterparts; for example, as 
far as we could check in the Literature, their simulation study 
only started at the end of the 1990'si^i. 

A few years before the existence of a BKT transition was 
proven for the LG-PR modeU 4 -, the same model had been in- 
vestigated by Monte Carlo simulation in the Grand-Canonical 
ensemble, and in the absence of a purely positional interac- 
tion, i.e. A = Oi^. Simulations were carried out for ft — 0.1 
and yu = -0.2, and showed that the BKT transition survives, 
even for mildly negative p, and that the transition temperature 
is an increasing function of p, in broad qualitative agreement 
with previous Renormalization Group (RG) studie a 17 ' 18 ' 19 . 

The Hamiltonian (Eq. ©) can be interpreted as describing 
a two-component system consisting of interconverting "real" 
(vk = 1) and "ghost", "virtual" or ideal-gas particles (v* = 0); 
both kinds of particles have the same kinetic energy, p denotes 
the excess chemical potential of "real" particles over "ideal" 
ones, and the total number of particles equals the number of 
available lattice sites (semi-Grand-Canonical interpretation). 
The semi-Grand-Canonical interpretation was also used in 
early studies of the two-dimensional planar rotator, carried 
out by the Migdal-Kadanoff RG techniques, and aiming 
at two-dimensional mixtures (films) of 3 He and 4 H o 17,18 ' 19 , 
where non-magnetic impurities correspond to 3 He. In addi- 
tion to phase separations, the results in Refsi 17,18,19 show that, 
in a regime of low fraction of non-magnetic impurities, the 
paramagnetic phase undergoes a BKT transition. 

Notice also that the above Hamiltonian (Eq. ©) describes 
a situation of annealed dilution; on the other hand, models in 
the presence of quenched dilution, and hence the effect of dis- 
order on the BKT transition, have been investigated using the 
PR mode l 20,21,22 ] 23 ! 24,25 and very recently the XY model^ 5 .; it 
was found that a sufficiently weak disorder does not destroy 
the transition, which survives up to a concentration of vacan- 
cies close to the percolation threshold. 

In this paper, we present an extensive study of the ferromag- 
netic LG-PR and LG-XY models, whose Hamiltonian can be 
explicitly written as 



A 



U<k) 



(uj,\Uk,i + Uj t2 u kt2 )] - fiN, (5) 



in order to gain insights into their critical behaviour. The mod- 
els are further simplified by choosing A = 0, i.e. no pure po- 
sitional interactions; Let us also remark that two-component 
spins are involved in the PR case, whereas XY involves three- 
component spins but only two of their components are in- 
volved in the interaction: in this sense the two models en- 
tail different anchorings with respect to the horizontal plane 



in spin space; moreover, an even greater variety of anchorings 
can be realized via the recently introduced generalized XY 
models 2 ^ 2 !. 

The rest of the paper is organized as follows: in Section 
|II]we discuss details our simulation procedure. Section HITl is 
devoted to the discussion of simulation results: we have found 
evidence pointing to a BKT transition, and used the relevant 
finite-size scaling theory to locate it. Possible effects of the 
chemical potential on the nature of the transition are discussed 
in Section[IV] which summarizes our results. 



II. MONTE CARLO SIMULATIONS 

A detailed treatment of Grand-Canonical simulations can 
be found in or via Refsi 15,28 ' 29 ; the method outlined here 
has already been used in our previous studies of other LG 
model a 16,28,30 . To avoid surface effects simulations were car- 
ried out on periodically repeated samples, consisting of V = 
Lr sites, L = 40, 80, 120, 160; calculations were carried out in 
cascade, in order of increasing reduced temperature T. 

The two basic MC steps used here were Canonical and 
semi-Grand-Canonical attempts; in addition two other fea- 
tures were implemente d 31 ' 32 : (i) when a lattice sites was vis- 
ited, Canonical or semi-Grand-Canonical steps were ran- 
domly chosen with probabilities !P can and Pgc, respectively; 
we used l P cdn /'PGC = n — 1, since spin orientation is de- 
fined by (« - 1) angles, versus one occupation number and 
(ii) sublattice sweeps (checkerboard decomposition) 3 ^ 3 ^; thus 
each sweep (or cycle) consisted of 2V attempts, first V at- 
tempts where the lattice sites was chosen randomly, then V/2 
attempts on lattice sites of odd parity, and finally V/2 attempts 
on lattice sites on even parity. Equilibration runs took be- 
tween 25 000 and 200 000 cycles, and production runs took 
between 250 000 and 1 000 000; macrostep averages for eval- 
uating statistical errors were taken over 1 000 cycles. Different 
random-number generators were used, as discussed in ReL^ 2 .. 

Computed observables included mean Hamiltonian per site 
and its temperature derivative (specific heat at constant p and 
V), density and its derivatives with respect to temperature and 
chemical potential, defined by 



H* = - (A) . 



and by the fluctuation formulae (see e.g. 29 ) 



Pt 



dp 
3T 



/--»' 



1 

VT- 



[<jVA>-<jV)<A>], 



dp 

p ^Tp 



C^v = 1 lfl(A) 
k B k B V 8T 



(6) 



(7) 



(8) 



(9) 



VT 1 



[(A 2 )-<A) 2 ]. (10) 
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We also calculated mean in-plane magnetic moment per site 
and in-plane susceptibility defined by 

M=i(VF?) (11) 

and 

*=JL<F-F>, (12) 

where 

F = 2^ v k (Kjfc,iei + M t> 2e 2 ) , (13) 
k 

taking into account only the in-plane components of the vec- 
tor spin. 

A square sample of V sites contains 2V distinct nearest- 
neighbouring pairs of lattice sites; we worked out pair oc- 
cupation probabilities, i.e. the mean fractions Rjk of pairs 
being both empty (R ee = /(l - vj)(l - v^), both occupied 
(Roo = (vjVk)), or consisting of an empty and an occupied 
site (R eo = (( 1 - Vj)Vk + (1 - v J t)v J ^). It should be noted that 

Ree + Roo + Reo = 1 ■ 

Short- and long-range positional correlations, were com- 
pared by means of the excess quantities 

R' 00 =\n{^,R': =R 00 -p\ (14) 

collectively denoted by R* (notice that these two definitions 
entail comparable numerical values). 

Quantities such as p, pj, p M and the above pair correlations 
R jk or R* oa can be defined as "fluid-like", in the sense that they 
all go over the trivial constants in the SL limit. Let us also re- 
mark that some of the above definitions (e.g. C„ y and p T ) 
involve the total potential energy both in the stochastic vari- 
able and in the probability measure ("explicit" dependence), 
whereas some other definitions, e.g. p M or the quantities R jk, 
involve the total potential energy only in the probability mea- 
sure ("implicit" dependence). 

III. SIMULATION RESULTS 

We start by discussing the outcome of the simulations per- 
formed for the LG-XY model. Results for a number of ob- 
servables, such as the mean energy per site H*, and density 
p (not reported here) were found to evolve with temperature 
in a smooth way, and to be independent of sample sizes; the 
derivatives C^y and pj (both plotted on FIG. Q3 showed a 
rather smooth trend with temperature, and a "sharp" peak at 
T ~ 0.6, around which the sample-size dependence of results 
became slightly more pronounced; notice that neither the spe- 
cific heat nor pj are exhibiting a divergence as a function of 
the temperature. 

At all investigated temperatures simulation results for M 
(not reported here) exhibited the expected power-law decay 




FIG. 1: Simulation estimates for the specific heat per site and for 
—pr, obtained with different sample sizes: circles: L = 40; squares: 
L = 80; triangles: L = 120; diamonds: 160. The value /i = 
was used in the present calculations. Statistical errors for C„ y, not 
reported, range between 1 and 3%; otherwise, here and in the follow- 
ing Figures, statistical errors fall within symbol sizes. 

with increasing sample size, i.e. they were well fitted by the 
relation 

lnM = -bilnL + bo, bi > 0, (15) 

where the ratio b\{T)IT was found to increase with tempera- 
ture; a spin-wave analysis of the SL-PR (Refsi22i2i and those 
quoted therein) predicts a similar sample-size dependence, in 
agreement with simulation results for the present model; fur- 
thermore, in the low-temperature limit b\ becomes propor- 
tional to T. 

The behaviours of specific heat and susceptibility (FIG. |2| 
suggest that the LG-XY model undergoes a BKT-like phase 
transition. According to the BKT theory, in the thermody- 
namic limit, the in-plane susceptibility, x, diverges exponen- 
tially while approaching the transition temperature in the 
high temperature region^ i.e. 

X~a x exp[b x (T -&)-*]• T -> © + (16) 

and remains infinite in the low temperature region T < 0. 

For a finite sample with volume V = Lr, the situation is dif- 
ferent, since the susceptibility obeys the constraint ^ < V/T; 
thus x is always finite, and its exponential divergence in the 
transition region is rounded, because of the important finite- 
size effects. Actually, results for ln^ versus temperature 
(FIG. |2]i were found to be independent of sample size when 
T > 0.62, and showed a recognizable increase with it (a 
power-law dependence of x on L) when T < 0.6. 

In the critical region, i.e. T ~ 0, the correlation length is 
of the same order as the linear size, L, of the system. In this 
case, the exponential divergence disappears in the critical re- 
gion, but a reminiscence of the divergence can still be found 
in the behaviour of x at higher temperatures, where the cor- 
relation length is still small compared to large linear system 
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FIG. 2: Simulation estimates for the logarithm of the magnetic 
susceptibility x versus temperature, obtained with different sample 
sizes; same meaning of symbols as in Fig. Q] The solid curve is ob- 
tained by fitting to equation d!6t with data corresponding to L = 160. 

sizes and finite-size effects can be neglected. We first fitted 
our MC results obtained for the largest sample size and for 
temperatures in the range 0.60 < T < 0.62 to eq. ( fT6l ). and 
obtained evidence of the exponential divergence of the suscep- 
tibility (see FIG. O, and an estimated transition temperature 
= 0.57 ±0.01, with a x = 0.123+0.005 and b x = 1.55 + 0.01. 

We checked our results by using data on a wider range of 
temperatures extending up to 0.70, and which yielded a con- 
sistent result. In the following we refine our results using a 
more elaborate method, namely the finite-size scaling theory. 
Details on this procedure can be found in or via referenced. 

In the temperature region where the singularity of \ is 
rounded the finite-size scaling theory holds, i.e. the corre- 
lation length £ is proportional to the size of the sample, £ ~ L. 
From the behaviour of the magnetic susceptibility x ~ ^ 
with r\ — 1/4-, we end up with 

X(®) ~ L 2 -\ (17) 

where we have omitted the corrections arising from the pres- 
ence of the background (analytic) contributions to the finite- 
size scaling. 

The temperature region around T = 0.57, was analyzed in 
greater detail, at first by carrying out a linear fit of ln^ versus 
InL and extracting rj from the slope; the values were found 
to be j](T) = 0.223 + 0.008, 0.244 ± 0.012, 0.255 + 0.007, 
0.273 + 0.005, 0.308 + 0.015 for T = 0.5700, 0.5725, 0.5750, 
0.5775, 0.5800, respectively. A non-linear square fit, based 
on equation ( fTTT i. was attempted as well, and yielded results 
consistent with these ones; it also proved convenient for result 
visualization to plot lnfjfZT 7 ' 4 ) versus InL (FIG. [3]). Thus the 
transition temperature is estimated to be © = 0.574 + 0.003, in 
agreement with the above mentioned result obtained by fitting 
the data of the susceptibility in the high temperature region; 
at this temperature we expect the value of 77 to be j to within 
statistical errors. This result shows that the presence of "va- 
cancies" in the sample reduce the transition temperature by 



TABLE I: Comparison of the transition temperatures, & PR , of the 
LG-PR model obtained in RefiAS and those obtained using the finite- 
size scaling presented in the present paper. 





0™ (Ref^) 


&p R (present work) 


-0.2 


0.73 + 0.01 


0.71+0.01 


0.1 


0.79 ± 0.01 


0.75 ± 0.01 



approximately 18% compared to the SL transition tempera- 
ture, but does not change its nature. 
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FIG. 3: Plot of In (^L- 7 ' 4 ) versus InL. Here the curve closest to a 
horizontal straight line signals the transition. 

The PR counterpart of the present model had been investi- 
gated in Ref. 15 , where a maximum sample size corresponding 
to L — 80 was used. Moreover, in that paper the transition 
temperature was estimated as the temperature where x started 
exhibiting a power-law dependency on L. Calculations had 
been carried out with two different values of the chemical po- 
tential, fi = -0.2 and fi = +0.1, respectively (see e.g. Table 
[D; the transition had been found to survive in both cases (even 
with a slightly negative value), and the transition temperature 
had been found to increase with increasing fi, in qualitative 
agreement with previous RG results^. 

For further comparison we reexamined the results obtained 
in the case of the LG-PR model in two ways. On the one hand, 
the above finite-size analysis was applied to the simulation 
data produced in Ref. - . The transition temperatures, reported 
in Table UJ corresponding to the two values of the chemical 
potential used were found to be lower than the ones reported 
previously. The behaviour of the ju-dependence of the tem- 
perature was similar for both methods. On the other, in order 
to achieve a better comparison with the present case, simula- 
tions were run anew for the case fi — 0, using same sample 
sizes as for the LG-XY counterpart (L = 40,80, 120, 160). 
The thermodynamic quantities of interest were found to be- 
have qualitatively in a similar fashion as those obtained in the 
framework of the XY model. 

A detailed investigation of the susceptibility via the above 
finite-size scaling procedure resulted in a BKT-like transition 
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with a transition temperature estimate of 0.733 ± 0.003, corre- 
sponding to a particle density of 0.924 ± 0.003. The transition 
temperature in this case was found to be 18% lower than the 
case of the planar rotator SL counterpart. 

As for other fluid-like quantities obtained within the frame- 
work of LG-XY, the plot of the derivative p^ (FIG. @]i against 
the temperature exhibited a smooth evolution, a broad max- 
imum at T ~ 0.625 and a comparatively weak sample-size 
dependence; pair occupation probabilities R JK were found to 
be essentially independent of sample size, similarly to quanti- 
ties such as H* and p; their results for the largest sample size 
L = 160 are presented in FIG. [4] as well; the three quanti- 
ties evolve with temperature in a gradual and monotonic way, 
and their plots suggest inflection points between T = 0.6 and 
T = 0.625, roughly corresponding to the maximum in p^. 
Short- and long-range positional correlations have been com- 
pared via the excess quantities R* og , whose simulation results 
for the largest sample size are reported in FIG. [5] showing a 
broad maximum well above the transition temperature; notice 
also that R* 00 remains rather small, reflecting the absence of a 
purely positional term in the interaction potential. In general 
it has been found that the obtained properties exhibit a recog- 
nizable similarity with with the counterparts for the isotropic 
LG-PR modet^. Some other remarks are appropriate, con- 
cerning the temperature dependence of fluid-like properties; 
let us recall that our potential model contains no explicit posi- 
tional interaction, i.e A — in Eq. ©, so that positional prop- 
erties are essentially driven by the orientation-dependent pair 
potential. Similarities and differences in the above plots for 
C^v and pj on the one hand (FIG. [U, and p M or R jk (FIG. |4]i 
on the other hand, can then be correlated with the "implicit" or 
"explicit" dependencies of the named quantities on the poten- 
tial energy (see remarks at the end of the previous section); the 
pair interaction energy changes most rapidly at T » 0.6, and 
this is reflected by narrow peaks in C^y and pj\ conversely, 
quantities with "implicit dependence" on the pair interaction 
energy are less affected by the change, as shown by a broader 
peak at higher temperature T =s 0.625. 



IV. CONCLUDING REMARKS 




FIG. 4: Simulation estimates for the three pair occupation proba- 
bilities R j K obtained for a two-dimensional sample with linear size 
L = 160, along with simulation estimates for p n , obtained with dif- 
ferent sample sizes; same meaning of symbols as in FIG.Q] 
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FIG. 5: Simulation estimates for the quantities R* o0 ; discrete symbols 
have been used for simulation results obtained with L = 160, and 
have the following meanings: circles: R' oa \ squares: R" a . 



In this paper, we have investigated two magnetic lattice- 
gas models living on two-dimensional lattices, i.e. PR and 
XY, respectively. We have studied both lattice-gas models, 
in the absence of a pure positional interaction, using semi- 
Grand-Canonical Monte Carlo Simulation. In the case of zero 
chemical potential a number of thermodynamic quantities ob- 
servable, including some characteristic of fluid systems, were 
computed. In general it has been found that the obtained prop- 
erties for both models exhibit a recognizable qualitative simi- 
larity. 

Results for the susceptibility yield consistent evidence of 
the existence of a BKT transition, suggesting the estimates 
®xy = 0.574 ± 0.003 and & PR = 0.733 + 0.003 for the tran- 
sition temperatures, with corresponding densities at transition 
p XY = 0.918 ± 0.004 and p PR = 0.924 ± 0.003, respectively. 
These results are consistent with those reported in Refi^ 5 . ob- 



tained in the framework of models describing quenched dilu- 
tion; the ratio of the transition temperatures of the LG mod- 
els to the SL ones is approximately 0.82. The peak of the 
specific heat falls closer to &bkt than for the SL counter- 
part; this can be interpreted as reflecting the fact that, since 
A = 0, the strengthening of short-range orientational correla- 
tions also brings about a more rapid increase of density, and 
hence, in turn, of H*. 

As for the possible /i-dependence of the present results, let 
us first notice that, since (dp/dp)ry > 0, it is reasonable to 
expect that %, and hence &lg, are increasing functions of p. 
On the other hand, the ground-state energy for the SL coun- 
terpart is — Vq per site, where Vq = 2; thus the lattice is es- 
sentially saturated when T < 1 and p exceeds a few multiples 
of Vq (say p ranging between 5 and 10Vo); at the other ex- 
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treme, when is negative and of comparable magnitude, one 
expects an essentially empty lattice, where no BKTLT should 
survive. Moreover, a BKTLT may still survive if // is negative 
but not too large in magnitude; this conjecture is supported by 
RG studies for the LG-P R 17 ' 18 ' 19 , and was confirmed by the 
simulation results obtained in Ref.— ; it also agrees with the 
conclusions of reference a 22 ' 23 ' 24 ; 25 for two-dimensional mag- 
nets with quenched dilution; for large negative //, the named 
RG treatments predict phase separation, i.e. a first-order tran- 
sition between a BKT and a paramagnetic phase, and even- 
tually the disappearance of the BKT phase. We expect to in- 
vestigate the ju-dependence of the transition temperature in a 
more detailed way in future work. 
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